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ABSTRACT 



Betelgeuse, the bright, cool red supergiant in Orion, is moving supersonically relative to the local interstellar medium. The star emits 
a powerful stellar wind that collides with this medium, forming a cometary structure, a bow shock, pointing in the direction of motion. 
We present the first 3D hydrodynamic simulations of the formation and evolution of Betelgeuse's bow shock. The models include 
realistic low-temperature cooling and cover a range of plausible interstellar medium densities of 0.3 - 1.9 cm~^ and stellar velocities 
of 28 - 73 kms~^ We show that the flow dynamics and morphology of the bow shock difl'er substantially because of the growth of 
Rayleigh-Taylor or Kelvin-Helmholtz instabilities. The former dominate the models with slow stellar velocities resulting in a clumpy 
bow shock substructure, whereas the latter produce a smoother, more layered substructure in the fast models. If the mass in the bow 
shock shell is low, as seems to be implied by the AKARI luminosities (~3xlO~^ M©), then Betelgeuse's bow shock is very young and 
is unlikely to have reached a steady state. The circular nature of the bow shock shell is consistent with this conclusion. Thus, our 
results suggest that Betelgeuse only entered the red supergiant phase recently. 
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1. Introduction 

Massive stars have a tremendous impact on their circumstellar 
environment. Betelgeuse (a Orionis, HD 39801, HIP 27989), 
the prototype red supergiant (RSG) and the brightest M su- 
pergiant in the sky (M1-M2 la-Iab), is no exception. Although 
Betelgeuse's initial mass is uncertain, it is thought to be >1OM0 
(see Table. [T]). Thus, the star is likely in the core helium-burning 
phase and is approaching its final stages of evolution - ultimately 
exploding as a core-collapse supernova and enriching the inter- 
stellar medium (ISM) with radiation, mechanical energy, mass, 
and heavy elements. Owing to its high luminosity, -lO^L©, and 
powerful stellar wind, Mwind~3xlO"^ Moyr"^ (see Table. [B, the 
star is already injecting copious amounts of mass and energy into 
the ISM. 

For stationary stars, the wind-ISM interface is expected 
to be a spherical shell as shown, for example, by the ana- 



lytic an d hydrodynamic models of Weav er et al.l (1 1977) ) and 
iGarcia -Segura et al. (1996), respectively. Betelgeuse has an av- 
erage radial (heliocentric) velocity of Vrad = +20.7 ± 0.4 kms"^ 
away from the sun, and the tangential velocities are Vacoss 
= 23.7(7) rpc1/20Q ) kms"^ and Vs = 9A(D [pc]/200) kms"^ 
dHarper et al."2008', and references therein). The distance D to 
Betelgeuse is not yet well established, the most recent estimate 
being 197+45 pc (see lHarper et al.|[2008L for a detailed discus- 
sion). For this range of distances, Betelgeuse is moving super- 
sonically relative to the local ISM, and the collision of its stel- 
lar wind with this medium has resulted in the formation of a 
cometary structure, a bow shock, pointing in the direction of mo- 
tion (at a position angle of 69.0° east of north). Although it is the 
only known RSG runaway, theoretical models predict that up to 
30% of RSGs can be runaway stars (Eldridge et al. 201 1). 

Early observations of the circumstellar medium around 
Betelgeuse suggested there is extended emission in the far- 




Fig. 1. Bow shock and 'bar' of Betelgeuse. [Left] The 60 
Jim IRAS observation retrieved from the IR AS Galaxy Atla s 
(IGA) Image ServeiEI(for details about IGA, see lCao etal]ll997h . 
[Right] The high-resolution, composite AKARI image with 
65 yum, 90/ im, and 140 pim emission in blue, green, and red, re- 
spectively (lUeta et al.ll2008.) [Credit: AKARI MLUES team]. 



infrared (e.g. IStencel et al.lll988l:r^ng et al.|[T99 3bl). However, 
because the star is so bright at these wavelengths, it was only in 
1997 that the bow shock was imaged success fully using IRAS 
Sit 60 and 100 iim (iNoriega-Crespo et aDll997h . More recently, 
"Ueta et alj (l2008l) have imaged the bow shock with AKARI at 
higher spatial resolution in the 65, 90, 140, and 160 yum bands. 



Send offprint requests to: shazrene@saao.ac.za 
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The bow shock is a smooth, rather circular arc that has a thick- 
ness of ~1.5 arcmin and an inner radius of ~5 arcmin (see 
Fig. [B. Also detected by both IRAS and AKARI is a mysteri- 
ous 'bar' structure, of unknown origin, located just ahead of the 
bow shock arc. 

Many of the observed properties of the bow shock can be de- 
rived in terms of parameters of the interacting flows, making this 
an invaluable probe of the physical conditions in both the stel- 
lar wind and the ISM. The characteristic size of the bow shock 
shell, known as the stand-ofl' distance, Rso, is determined by the 
location where the ram pressures of the two 'flows' are in equi- 
librium: 

where pisM and pw are the density of the ISM and stellar wind, 
respectively; is the velocity of the star with respect to the ISM, 
and Vw is the stellar wind velocity. Assuming spherical mass loss 
from the star, 

PwVw = Mw/47^r^ (2) 

and rearranging the resulting equation yields the well known so- 
lution: 

\ 47rpisMv: 

Making the additional assumption of momentum conservation 
and assuming that the material mixes and cools instantaneously 
(so that the dense sh ell has negligible thickness, i.e. the thin- 
shell approximation) , Iwilkini |l996) derived an analytic expres- 
sion for the shape of the bow shock: 

R(0) = 7?so cosec ^J3(l-0cote) , (4) 

where 6 is the polar angle measured from the axis of symmetry. 

Utilising the above analytic models and current estimates 
of Betelgeuse's wind properties and distance (see Table. [B, 
lUeta et al.l (l2008l ) derived a space velocity of = 40 n^ 
kms"^ with respect to the local ISM. Estimates of the ISM 
density, nn, range from 0.3 cm"^ for a flow em anating from 
the Orion Nebula Complex (ONC) (iFr isch et a l. 1990) to 1.5 - 
1.9 cm~^ if the origin of the flow is the Orion OB 1 associa- 
tion (lUeta et al.ll2'008[[2009) . Given this range of ISM densities, 
Betelgeuse's space velocity with respect to the ambient medium 
is therefore likely to be between 73 kms"^ and 28 kms"^ If 
we assume strong shock conditions, the post-shock temperature 
corresponding to these stellar velocities, is at least 10000 K. To 
date, the bow shock has only been detected in the far-infrared, 
where the bulk of the emission is probably caused by dust grains 
with an uncertain contribution from oxygen and carbon fine- 
structure lines. 

There are several multi-dimensional hydrodynamic mod- 
els of wind-ISM interac tions with relevant parameters for 
Betelgeuse. For example, iGarcia-Segura et aP |l996) have in- 
vestigated the interaction of the ISM with a wind ejected from a 
stationary star evolving from the main sequence to the RSG and 
then Wolf-Rayet (WR) phase. They found that the RSG wind 
had a significant eff'ect on the subsequent WR phase. The same 
conclusion was reached by Brighenti & D'Ercole (1995) who 
studied the circumstellar medium (CSM) resulting from a mov- 
ing star in 2D. As expected, in this case a bow shock arc was 
formed rather than a spherical shell. The implications of an RSG 
phase for the WR progenitors of gamma-ray bursts was also 
discussed in a similar inv estigation by yan Ma rie et al. ( 2006). 
IWareing et ^ (l2007ah and lVillaver et al.r(l2Q03h calculated a se- 
ries of 3D and 2D hydrodynamic models, respectively, of the 



Table 1. Summary of Betelgeuse's parameters. 



Parameter 


Value 


Refs. 


Mass (M=,) 


11.6!^-^, 15_20 Mo 


[1,2] 


Temperature (T^) 


3 300K 


[2] 


Radius (R^) 


950 - 1 200 Ro 


[2] 


Luminosity (L*) 


0.9-1.5 xlO^Lo 


[2] 


Wind velocity (v^) 


17 ± 1 kms-i 


[3] 


Mass-loss rate (M^) 


2-4x10-6 Moyr-i 


[4] 


ONC ISM density (^h) 


0.3 cm"^ 


[5] 


ISM density (^h) 


1.0-1.9 cm"^ 


[6] 


Stellar velocity (v*)^ 


73 - 28kms-i 


[6] 


Stand-off' distance (Rso) 


(8.5 + 1.9)x 10^^ cm 


[6] 


Distance (D) 


197 ± 45 pc 


[7] 



Notes, a. Stellar velocity with respect to the local ambient ISM. 
References. (1) Neilson et al. (2011). (2) Smith et al. (2009) an d refer - 
ences therein. (3) Bernat et al. (1979). (4) Noriega-Crespo et all (119971) . 
(5) Observed ISM densit y of the Orion N ebulae Complex (ONC), 
Frischetal. (1990). (6) lUeta et all (l2008h . bow shock fitting. (7) 
Harper et al..(.2008.) . 

interaction of a stellar wind from an asymptotic giant branch 
(AGB) star with the ISM and the subsequent structure formed 
during the planetary nebula phase. Other models of the CSM 
around fast-moving AGB star s include 3p and axisymmetric 
models of Mira's bow shock ("Raga et al." "2008'; 'Esquivel et al.l 
2010, respectively). More recently, van Marie et al. ( 2011) have 
modelled Betelgeuse's bow shock with 2D high-resolution sim- 
ulations that include a simple dust tracking scheme. They show 
that the ffow for grains of various sizes diff'ered and can have 
important consequences for interpreting infrared observations of 
bow shocks. 

We build on the work above by including a realistic treat- 
ment of the thermal physics and chemistry in 3D models and 
consider a broader range of parameters. The numerical method 
and model parameters are described in Sect. O In Sect. O we 
present an adiabatic model that serves as a code test and a point 
of reference for the simulations that include radiative cooling. 
The ffow characteristics, morphology, and shell properties of the 
latter are described in Sect. |4l Finally, the implications of this 
work, particularly with regard to Betelgeuse's bow shock mor- 
phology, shell mass and age are discussed in Sect. \5\ 

2. Numerical method 

2.1. Smoothed particle hydrodynamics 

Smoothed particle hydrodynamics (SPH) is a Lagrangian 
method in which particles behave like discrete ffuid elements. 
However, each of their ffuid properties, e.g. density, temperature, 
or velocity, is the result of mutually overlapping summations and 
interpolations of the same properties of neighbouring particles. 
In this way, a continuous ffuid is realised; i.e. physical ffuid prop- 
erties that vary smoothly over all points in space can be defined 
from interpolations on a finite number of particles. 

SPH has been utilised to study a wide range of problems, not 
only in astrophysics, but also in other fields, e.g. engineering. 
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Table 2. Summary of the parameters used for the different com- 
puted CSM evolution models. 



Model 




hh 


TiSM 




Comment 




(kms-i) 


(cm"^) 


(K) 


(0.01 pc) 




A 


28.8 


1.9 


650 


0.3-2 


Cooling 


B 


32.4 


1.5 


1000 


0.2-2 


Cooling 


C 


39.7 


1.0 


1600 


0.2-2 


Cooling 


D 


72.5 


0.3 


8 000 


0.1-2 


Cooling 


Dl 


72.5 


0.3 


8 000 


0.5-4 


Cooling 


Cl650 


39.7 


1.0 


650 


0.6-3 


Cooling 


Cl1600 


39.7 


1.0 


1600 


0.6-3 


Cooling 


Clsooo 


39.7 


1.0 


8 000 


0.6-3 


Cooling 


CLf 


39.7 


1.0 


1600 


0.6-3 


/2~0.18 


Bad 


32.4 


1.5 


100 


1-5 


Adiabatic 


Ah 


28.8 


1.9 


650 


0.03 - 0.6 


Cooling 



Notes. 1 . The range h is an average of the smaller smoothing lengths in 

the bow shock to the typical ISM value. 

2. / is the fraction of molecular to atomic hydrogen. 



2.2.1. The ISM 

The ISM is assumed to be homogeneous. Although a cartesian 
grid of particles is the simplest approach to producing such a 
smooth medium, it also results in unwanted artefacts along pre- 
ferred grid directions. To overcome this problem, we use a 'glass 
tile', where particles are randomly distributed in a periodic box 
and evolved with the sign of gravit y reversed unt il an equilib- 
rium configuration is reached (see Springel"2005', for details). 
The relaxed glass tile is then scaled and cropped to the required 
dimensions for the simulation. 

We model a range of ISM densities, nn = 0.3, 1.0, 1.5, and 
1.9 cm"^, with corresponding stellar velocities = 73, 40, 32, 
and 29 km "\ respectively (see Table. [2l). These number densities 
lie at the boundary between typical values expected for either a 
warm or cold neutral ISM, so we assume tempe ratures based 
on the phase diagram of the standard model of IWolfire et al.l 
( Il995h (e.g. their Fig. 3d). The temperatures, Tjsm, are 8 000, 
1 600, 1 000, and 650 K, respectively. The mass of each particle 
is determined by the volume of the simulation domain, Vbox, the 
number of ISM particles, Msm, and the gas number density as 
follows: 

^particle ~ -kj ' \^) 

^ViSM 

where ji = 1.4 is the adopted mean molecular weight and mu 
the mass of atomic hydrogen. The chosen particle mass eff'ec- 
tively determines the stellar wind resolution, i.e. the number of 
particles required to reproduce the observed mass-loss rate, as 
described in the section below. 



Thus, the SPH equations are often adapted and optimised for 
a specific field or application. We briefly outline the particular 
SPH formulation implemented i n the code, wh i ch serves as a 
basis for this work, GADGET-2 (ISpringell[2Q05t ISpringel etaP 
12007), and refer the reader to detailed reviews of the method in 
,Monaghan (1992), Rosswog (2009), and Springel (2010). 

The momentum equation is derived based on an 'entropy for- 
mulation' for SPH, which conserves entropy, energy, and mo- 
mentum by construction ( Springel & Hernquist 2002). In this 
formulation, a function of the entropy, A(s) = Pjp^, where y 
is the ratio of specific heats, P the pressure, and p the density, is 
evolved instead of the internal energy, e. To model shocks, vis- 
cosity terms that transform kinetic motion of the gas into internal 
energy are also included in the momentum equation. When the 
fluid is compressed, the viscosity acts as an excess pressure as- 
signed to the particles in the equation of motion. 

GADGET-2 also includes s elf-gravity and utilis es a hierarchi- 
cal Barnes and Hut oct-tree (iBarnes & Hutlll986l) algorithm to 
calculate the gravitational accelerations. The simulations carried 
out in this work are characterised by a broad dynamic range in 
densities. To obtain better and more eflftcient spatial and time 
resolution, particles have adaptive smoothing lengths and are 
evolved with individual and adaptive timesteps. 

2.2. Model setup 

For numerical convenience, we choose a frame in which the star 
is stationary and located at the origin (x, y,z = 0, 0, 0) of a rect- 
angular box. The ISM flows past the star in the direction of the x 
axis, interacting with the stellar wind as it does so. Two difl'erent 
flags are utilised to distinguish between ISM and wind particles. 



2.2.2. The stellar wind 

Although the extended envelope of Betelgeuse is thought to be 
inhomogeneous and asymmetric close to the star, on much larger 
scales, and more importantly on the scales that we a re interested 
in, the outflow is largely sphericall y symm etric (ISmith et al.l 
120091) . Furthermore, as discussed in [Harped (l2010h . the mass- 
loss mechanism for M supergiant winds is not well understood. 
Given these uncertainties, we do not model the wind acceleration 
in detail, and instead assume that the wind is launched isotropi- 
cally and with constant velocity. 

The particle mass calculated above then determines the num- 
ber of particles, A^wind, that are injected in an interval A^^ind to 
produce a constant mass-loss rate of 3.1 x lO'^Moyr"^ The 
wind particles are injected as spherical shells of typically a few 
hundred particles "cut" from the glass tiles described above. 
Each shell is randomly rotated in the x,y, and z directions. 
Although this reduces the smoothness of the outflow, it pre- 
vents the occurrence of numerical artefacts. The particles are 
injected at a radius, 7?inner~10^^ cm, with a constant wind ve- 
locity, Vw~17kms"^. The initial temperature of the particles is 
^wind~l 000 K, the temperature at Thinner d erived from models of 
Betelgeuse's circumstellar envelope, e.g. iGlassgold & HugginsI 
(1986) and Rodgers & Glassgold (19911 

2.2.3. Cooling 

Radiative cooling plays an important role in determining the 
temperature structure and emission from the bow shock. Given 
the large unce rtainties in the proce sses involved, we adopt the 
approach of Smith & Rosen (2003*) rather than model the cool- 
ing in detail. Their cooling curve consists of analytic solutions 
and fits for various coolants based on d etailed calculations de- 
scribed extensively in the appendices of iSmith & Rosed (l2003h 
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Table 3. Summary of the combined cooling processes. 



Label 


Process^ 


^atomic 


Atomic cooling 


Aci,cii 


Carbon fine- structure cooling 


AcO(H) 


CO vibrations excited by atomic hydrogen 


^C0(H2) 


CO vibrations excited by molecular hydrogen 


AcO(r) 


CO rotational transitions collisionally 




excited by molecular and atomic hydrogen 


AcR 


Cosmic ray heating 


Agrain 


Dust grain cooling/heating 


AH2(diss) 


Dissociation of molecular hydrogen cooling 


AH2(ref) 


H2 reformation heating 


AH2(r-v) 


Molecular hydrogen rotational and 




vibrational cooling 


AH20(r) 


Water rotational cooling 


Ah20(H) 


Water vibrations excited by collisions 




with atomic hydrogen 


Ah20(H2) 


Water vibrations excited by collisions 




with molecular hydrogen 


Ao(63//m) 


Oxygen non-LTE cooling 


AoH 


OH cooling 



Notes . 1. The above processes are described in detail in lSmith & RosenI 
(l2003h and references therein. 



and references therein, e.g. rotational and vibrational transitions 
of H2, CO and H2O, H2 dissociative cooling and reformation 
heating, gas-grain cooling/heating, and an atomic cooling func- 
tion that includes non-equilibrium eff'ects (see Table. [3]). 

First the non-equilibrium molecular and atomic fractions 
of hydrogen are c alculated according to the prescription of 
ISuttner et^ll (1 19971) . and a limited equilibrium C and O chem- 
istry is used to calcul ate the CO, OH, and H2O abundances 
(ISmith & RosenI l2003l) . We adopt standard solar abundances 
of hydrogen and helium, and set the abundances of carbon 
and oxygen to a^c ^2 .5x1 0""^ and ;^o~6.3xlO"'^, respectively 
(iLambert et al.l 1 19841) . We assume that the initial fraction of 
molecular to atomic hydrogen, / = ^h2/^h, in the Betelgeuse 
outflow is small, / = 0.001 (on a scale where is atomic and 0.5 
is fully molecular). The hydrogen in the ISM is assumed to be 
atomic. The dust cooling function depends on the temperature of 
th e grains, and we se t this temperature to 40 K, the value derived 
bv lUeta et all (l2008h . The change in specific internal energy, 6, 
due to radiative cooling is then calculated semi-implicitly for 
each particle, /: 

,(..1) ^ ^ ^ ^J^Ll^At , (6) 

where A and T are the volumetric cooling and heating rates in 
ergscm"^ s"^ respectively; is the timestep; a nd n and n -\- I 
are the current and new timesteps, respectively (ISpringel et al.l 



simulations of Betelgeuse' s bow shock 









TCA/r 




RSG 











20 40 60 80 100 

Time [1000 years] 

Fig. 2. Position of the shock along the x axis (defined by the 
point with maximum temperature) as a function of time for the 
ISM (dashed line) and the stellar wind (dotted line). The black 
solid line demonstrates the trajectory for a Vg = 10 km s"^ shock 
velocity. All velocities are measured with respect to the station- 
ary frame of the star. 



I2001h . The change in specific internal energy is then converted 
to a change in the entropy of each particle. 

We do not include a treatment of dust formation and de- 
struction, depletion of elements and molecules on to grains, 
the eff'ects of magnetic or UV background radiation fields, den- 
sity inhomogeneities in the stellar wind, and ISM. Collisional 
ionisation, likely important at high temperatures is not explic- 
itly included; instead, hot gas, >10'^ K, cools primarily by 
atomic line emission ( Sutherland & Dopital 1 19931) and thermal 
bremsstrahlung. Although these processes are likely to have an 
eff'ect on the bow shock properties, both in terms of radiative 
transfer and cooling, the computational cost of their inclusion in 
3D models is prohibitive. 

In the following sections, we explore and discuss possible 
solutions for Betelgeuse's bow shock with 3D numerical mod- 
els (see Table. O. We simulate four basic models, A-D, with 
medium resolution, and with ISM number densities, nn = 0.3, 
1.0, 1.5, and 1.9 cm"^, and stellar velocities = 73, 40, 32, and 
29 km"\ respectively. Several lower resolution models are also 
calculated in order to study the long-term evolution of model D 
(model Dl) and the eff'ect of varying the assumed ISM tempera- 
ture (models Cl65o, Cli600, Clsooo) and the initial molecular hy- 
drogen fraction (model CLf). Model Ah has the same parameters 
as model A, but is computed with higher resolution. Finally, we 
have also computed an adiabatic model (model Bad), discussed 
in the next section, which provides a useful basis both for com- 
parison to previous studies and to our models with radiative cool- 
ing (Sect.il). 



3. Numerical test: Adiabatic model 

We tested the numerical set-up with model Bad, a purely adi- 
abatic model, i.e. no heat sources or sinks, and y = 5/3. The 
simulation begins at the start of the red giant phase, with the 
star losing 3.1 x 10"^Moyr"^ Moving at ~32kms"\ the star 
traverses a cold ISM with Tism = 100 K and density nu = 1.5 
cm ^. 
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X [pc] X [pc] 

Fig. 3. Adiabatic model, Bad, of the interaction of Betelgeuse's dense stellar wind with the ISM after t ^ 100000 years. We plot 
cross-sectional profiles along the symmetry axis (the xy plane) for p, the gas density [top, left]; T, the temperature [top, right]; the 
velocity vector field plotted over |v|, the gas velocity [bottom, left]; and the smoothing length, h [bottom, right]. Approximately 
1 X 10^ ISM particles filled the simulation domain and -6x10^ particles were injected to produce the isotropic stellar wind. 



The stellar wind collides with the ISM, and material accu- 
mulates at the contact discontinuity (the surface separating the 
two flows), where part of the kinetic energy of the gas is ther- 
malised. The heated ISM and wind material expand outwards 
from either side of this surface, the former moves into the ISM 
and is called the forward shock, and the latter into the stellar 
wind and is known as the reverse shock. The expansion of these 
two regions is shown in Fig.O The position of the reverse shock 
(defined as the point of maximum temperature for the RSG wind 
particles along the x axis) changes slowly and reaches its average 
position, — 0.25 pc, after 50000 years. The forward shock (de- 
fined as the point of maximum temperature for the ISM particles 
along the x axis) expands at 10 km s"^ during that time, and then 
decelerates to a position of -0.48 pc. The bow shock reaches a 
steady state after 80 000 years, and has an average width of -0.2 
pc. 

The density, temperature, and velocity in the bow shock are 
shown in the cross-section in Fig. [3j After 100000 years, the 
maximum width of the tail is approximately 3 pc (see Fig. [3] 
[top, right]). The average temperature in the tail is approximately 
4000 K; however, the cylindrical core of the tail is filled with 
much hotter, -15 000 K, gas spanning -0.8 pc in diameter. This 
hot gas is advected downstream from the forward shock, and 
flows towards the central part of the tail, filling the low-pressure 
void (beyond the RSG stellar wind) evacuated by the moving 
star (Fig. [3] [bottom, left]). 

While the cometary tail is largely composed of ISM mate- 
rial at early times, both the RSG wind and ISM contribute sim- 
ilar amounts of mass to the bow shock head. Defining material 
in the bow shock shell as shocked regions with T > 2 000 K 
and density above that in the ambient medium, we see in Fig.|4l 



:3 




20 40 60 
Time [1000 years] 



Fig. 4. Growth of the mass in the bow shock shell with time. 
The shell is defined as shocked regions with T > 2 000 K and 
density above that in the ambient medium. The lines show the 
total shell mass (black), the contributions from the ISM (blue) 
and the RSG wind (red) for x < Rso or 6 < 115°, whereas the 
points represent the same parameters except for shell material 
withx<0or6><90°. 



that after -40 000 years, the RSG wind contributes slightly more 
mass to the bow shock head (6 < 90°, i.e. x < 0) than the ISM, 
and it takes -90 000 years for this to occur for ^ < 1 15°, i.e. for 
X < Rso- The mass accumulates in the bow shock almost linearly 
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Fig. 5. Particle density [top] and temperature [bottom] profiles 
for model Bad in a thin slice, -0.03 <z< 0.03 and -0.03 <y < 
0.03, along the x axis. ISM and RSG particles are blue and red, 
respectively. 



after -30000 years and only reaches a steady state for 6 < 90° 
after ~80000 years. An estimate for a typical shell mass can be 
obtained from 

M^Rso 



At: 



shell 



(7) 



which yields 0.05 M©. This expression assumes the mass is dis- 
tributed over 471 sterad in a spherical shell and only takes the con- 
tribution from the stellar wind into account. Thus for < 90°, the 
expected RSG shell mass would be 0.025 M©, i.e. half of Alsheii- 
With a significant contribution from the wings of the paraboloid 
(the shell radius is greater than Rso for all 6> > 0°), the RSG shell 
mass obtained from the simulations is twice the theoretical value 
(red points in Fig.|4l). Including the mass from the ISM, the total 
mass in the bow shock head after 100000 years is ^0.1 M©. 

3.1. Numerical resolution and smoothing 

The smoothing length, h, determines the scale on which the 
above fluid properties are smoothed, so it efl'ectively defines the 
resolution of a simulation. In GADGET-2 the smoothing length of 
each particle is allowed to vary in both time and space, thus the 
full extent of SPH's powerful Lagrangian adaptivity is achieved. 
As shown in Fig. [3] [bottom, right], the resolution in high- 
density regions (e.g. in the bow shock) is naturally increased 
with large particle numbers, hence smaller smoothing lengths 
(/z<0.01pc), and little computational time is wasted evolving 
voids, which can be described by fewer particles with longer 
smoothing lengths (/z~0.08pc). 




X [pc] 

Fig. 6. Particles in the xy plane with the ISM in blue, and the 
stellar wind in red. The analytic solution (Eq. (4]) is shown in 
black. 



Unfortunately, this smoothing process also means that the 
shock that results from the collision of the ISM and the stel- 
lar wind is broadened over a few particle smoothing lengths. 
To minimise this broadening, it is essential that simulations 
employ as large a number of particles as possible to achieve 
shorter smoothing lengths and thus sharper discontinuities. In 
this simulation, we model the ISM with ~10^ particles and the 
stellar wind with ~10^ particles, achieving a typical smoothing 
length of 0.03 pc or better in the regions of interest. More im- 
portant, however, despite the shock broadening, the post-shock 
fluid properties are still described well; e.g. as expected from the 
Rankine-Hugoniot jump conditions, the density jump is a fac- 
tor of 4 (see Fig. [5] [top]). Assuming that all the kinetic energy 
in these winds is thermalised, we can derive an estimate for the 
expected temperature of the shocked gas. 



_ 1 jimnv^ 
'~2 



(8) 



where is the Boltzmann constant. For the parameters assumed 
here, this equation yields 19 000 K for the shocked stellar wind 
and 92 000 K for the shocked ISM temperature, so very similar 
to the results obtained in the simulation (see Fig. [5] [bottom]). 

3.2. Development of instabilities 

The shear produced by the relative motion of the shocked ISM 
and shocked stellar wind regions (as shown in the vector plot 
in Fig. [3] [bottom, left]), excites Kelvin-Helmholtz (K-H) insta- 
bilities at the con tact discontinuity. Utilising the prescriptions 
described by Brig henti & D'Ercold (Il995h (Eqs. on p55), we 
find that the K-H growth time scale is at least an order of mag- 
nitude smaller than the one required for the Rayleigh-Taylor 
(R-T) instability. Thus, it is the K-H rolls that are advected 
downstream, and the characteristic R-T 'fingers' that appear in 
Brighenti & D'Ercole (1995) (their Fig. 2 [top panel]) are ab- 
sent from our simulation. The main reason for the diff'erent be- 
haviour is that they utilise an order of magnitude higher mass- 
loss rate and a much slower wind velocity, both of which lead 
to a much stronger density contrast between the ISM and wind, 
which favours shorter growth scales for R-T instabilities. 

The development of K-H in stabilities in our m odel may be 
surprising. Previous studies, e.g. lAgertz et al.l (120071) . have found 
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that instabilities are suppressed in 'standard' SPH, owing to dif- 
ficulty estimating gradients in regions with strong density con- 
trasts. Several methods have been proposed to address this issue. 
One approach is to include artificial conductivity, which acts as 
an additional pressure across the contact d iscontinuity, facilitat- 
ing m ixing between the two layers (|Pricdl2QQ8l: IWadsley et al.l 
l2008h . However, at present the formulation is not consistent 
with including se lf-gravity and tends to be highly dissipative 
(IValcke et al.l 120101) . The other app roach, which was the result 
of a detailed study by IValcke et al.l (2010), is to minimise parti- 
cle disorder, which prevents the 'oily' problem (the development 
of an artificial gap between the fluids) and allows the gas to mix. 
They achieve this by using a higher order smoothing kernel to 
prevent particle clumping, but this kernel does not conserve en- 
ergy as well as the cubic spline kernel. Given the caveats to the 
above solutions, we opted to use glass initial conditions instead, 
i.e. to start from a very smooth particle distribution with mini- 
mal particle disorder. Thus, strong K-H instabilities do develop, 
despite the factor of 10 in the density contrast. 

3.3. Comparison to analytic models 

Although the contact discontinuity is distorted by the K-H insta- 
bility, the overall shape of this surface is in fairly good agreement 
with the one predicted by the lWilkinI (0^996) solution (Eq.[T]) for 
small polar angles (6), and the stand-oflf distance matches the 



position of the contact discontinuity well (see Fig. [6] [top]). In 
the thin-shell limit, the ISM and wind material are assumed to 
cool instantaneously and mix fully. In contrast, the forward and 
reverse shock of the adiabatic model have finite width, -0.2 pc, 
and as expected, there is little mixing between these two layers 
of gas. Overall, the adiabatic model demonstrated that the set-up 
and code can achieve results that are consistent with both theo- 
retical expectations and previous studies. 

4. Models with cooling 

The models that include radiative cooling exhibit a global struc- 
ture that is similar to the structure of the adiabatic model. In 
Fig. [71 we show the density, temperature, and emissivity for 
models A-D. The four regions described above can easily be dis- 
tinguished: the unshocked ISM and similarly unshocked, free- 
flowing stellar wind, separated by a layer of hot, shocked ISM 
(the forward shock) and shocked stellar wind (the reverse shock). 
Based on their flow characteristics and morphology, models A- 
D can be grouped into two classes: the 'slow' models, A-C, with 
ISM densities >1 cm"^ and peculiar velocities <40kms"^ and 
the 'fast' model, D, with an ISM density of 0.3 cm"^ and a stellar 
velocity ~73 km s"^ 

The fluid properties of the slow models show significant de- 
partures from the adiabatic case (see Fig. (8]); the most obvious 
being that the radiative cooling results in lower post-shock tem- 
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Fig. 8. Particle density [top], temperature [middle], and velocity 
[bottom] profiles for models B [left] and D [right] in a thin slice, 
-0.02 < z< 0.02 and -0.02 <y < 0.02, along the x axis. In all 
plots, the ISM and RSG particles are blue and red, respectively. 



ISM 
RSG 

V = 10 




Model B 



20 40 60 

Time [1000 years] 

Fig. 9. Same as Fig.[2]except for models B and Dl. 



peratures, e.g. for model B is an order of magnitude lower 
than in model Bad for both the shocked ISM and RSG wind 
(see Fig. [8] [top, middle]). Furthermore, by reducing the thermal 
pressure of the gas, the strong cooling also enables further com- 
pression in the shock, resulting in greater post-shock densities 
(see Fig. [8] [top, left]) than in model Bad- Consequently, the time 
scale for the growth of R-T instabilities is reduced and, similar 






log TA/Xdvn 

Fig. 10. Ratio of the cooling and dynamical time scales for mod- 
els Ah [left] and D [right]. If the cooling time scale is much 
greater than the dynamical time scale, i.e. ta » Tdyn, the gas 
behaves adiabatically, and if ta « Tdyn, the gas is approxi- 
mately isothermal. 



to the models of iBrighenti & D'Ercold (Il995h (their Fig. 2.), R- 
T 'fingers' develop faster than the K-H rolls that characterised 
the adiabatic case (see Secs. l3.2l and l4.1.ll) . Owing to the lower 
thermal pressure and mixing of the ISM and RSG gas, both the 
bow shock and the tail are narrower in the slow models; e.g., the 
model B width is approximately half that of model Bad- 

The forward shock of the fast model shows similar charac- 
teristics to model ^ad; the combination of a large stellar velocity 
through the ISM, and an initially hot and lower density ISM, 
means that the fast-moving gas does not have enough time to 
cool beyond its initial condition. The fast model achieves the 
same shocked ISM temperature, Tism, as the adiabatic case in 
which the star was moving approximately half as fast (see Fig. [8] 
[middle panel]). Furthermore, this hot gas also flows along the 
bow shock and eventually settles in the core of the tail along the 
X axis, as it did in the adiabatic model. However, the shocked 
RSG wind temperature is approximately 1 000 K, much lower 
than the adiabatic counterpart and similar to the value for the 
slow models. 

The evolution of the forward and reverse shock positions 
for models B and Dl, i.e. the position of maximum tempera- 
ture along the symmetry axis, is shown in Fig. [9l Not only are 
the oscillations in the shock position less violent than the adia- 
batic case, an equilibrium shock position is attained more rapidly 
for the models with radiative cooling, after only 40000 years. 
The strong cooling in both the forward and reverse shocks in 
the slow models results in the excitation of R-T type instabilities 
(see below), which causes the gas to mix, hence the small dif- 
ference between the positions of the RSG and ISM components. 
In contrast, there is little mixing between the two layers in the 
fast model's bow shock, the forward shock has significant width, 
thus the hottest gas extends much further into the impinging ISM 
flow. 

The typical cooling time scale in the bow shock for x < 
is -10000 years, and is of the order of both the characteris- 
tic dynamical time scale for the shocked ISM, 7?so/v*, which is 
-9400, 8 300, 6 800, 3 700 years for models A-D, respectively, 
and the characteristic time scale for the shocked wind, 7?so/Vw, 
which is approximately 17 200 years. In Fig.[TOlwe compare the 
radiative cooling and characteristic dynamical time scales of the 
gas: 



Ta = (A/pe) 



-1 



(9) 
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Time [1000 years] 





t = 73000 yrs^ 











Fig. 11. [Top] Ratio of 7?(0°)/7?(90°) as a function of time for 
models B (magenta) and D (blue) compared to the analytic value 
in black. [Bottom] Hydrogen column density (on a logarithmic 
scale) for model B after 43 000 years [left] and 73 000 years 
[right]. The analytic solution for the shape of the bow shock, 
given by Eq.|4l is shown in blackQ 



and 



Tdyn - ^SO/V , 



(10) 



respectively. For regions with ta « Tdyn, e.g. the bow shock 
arc where the gas is strongly decelerated, the gas cools strongly 
and behaves almost isothermally. The gas in the tail has a much 
longer cooling time scale, -30000 years (much greater than the 
cooling time scale for the shocked ISM) resulting in a largely 
adiabatic flow there. 



4.1. Bow shock morphology 

The shape of the bow shock changes with time: it is initially 
circular and becomes increasingly 'parabolic' until it reaches a 
steady state, at which point the glo bal morphology is described 
by Eq. m the analytic Twilkinl (1199 6) solution. The ratio of the 
shock position at angles 6> = 0° and ^ = 90°, rs(0°)/rs(90°), 
as a function of time, is shown in Fig. [TT] [top]. It takes sev- 
eral dynamical time scales for the bow shock to achieve the 
steady state value 7?(0°)/7?(90°)=l/ V3. The wind takes at least 
t = R(0)/v^ to reach the thin-shell position, ignoring the drag 
due to the moving environment, thus it takes -16000 years for 

6> = 0° (7?(0°)=7?so), while for 6 = 90° (7?(90°) = V37?so) it 
takes -28 000 years. In reality, the above time scales are longer. 
The drag is maximal for ^ = 0°, so it takes about 30000 years 
to reach R^o (e.g. Fig. O, and coincidentally it takes about the 



^ Movies demonstrating the evolving bow shock shape are included 
in the electronic version. 



same time to reach the solution position for = 90° since the 
distance traversed by the wind is greater, but the drag is reduced. 

It takes even longer for a steady state to be achieved for 
larger angles, 6 > 90°. In Fig. [TT] [bottom] we compare the thin- 
shell solution with model B at 43 000 years and at 73 000 years. 
Although the simulations coincide with the analytic solution in 
the latter, as expected, the agreement for large angles of 6 at ear- 
lier times is not good, and the same efl'ect is found for the fast 
model. The shape of the bow shock tail also changes with time, 
and younger systems show strong curvature in the bow shock 
tail. This is the case for both the slow and fast models; however, 
for the latter the curvature is quickly left far downstream due to 
the large space motion. 

The inclination angle, the angle between the apex of the bow 
shock and the plane of the sky, also has a marked efl'ect on the 
shape of the bow shock. In Fig. [12] we show the column densities 
for both the fast and slow models at difl'erent inclination angles. 
The bow shock becomes increasingly circular and the curvature 
in the tail also increases with larger inclination angles. At the 
same time, increasing the inclination angle reduces the density 
contrast between the peak at the apex of the bow shock and the 
rest of the cometary structure, making the detection and identifi- 
cation of highly inclined systems more difficult. 

4.1.1. Substructure: Instabilities 

The bow shock is unstable for all models, however, the slow and 
fast models clearly exhibit very different bow shock substruc- 
ture. The slow models show a thin, rather smooth outer shock 
and a contact discontinuity that is highly distorted (on a scale 
of -0.1 pc) by prominent R-T 'fingers' of the RSG wind (see 
Figs.[8l [T2l[T^IB.21 and App.|A|. The protrusions are even more 
developed in the high-resolution model Ah (see Figs. [13] and 
lA.ll) . In the column density plots, the small-scale R-T instabil- 
ities result in a clumpy/knot-like structure that becomes one of 
the dominant features of the bow shock at large inclination an- 
gles (see Fig. [12]). 

In contrast, the bow shock in the fast model is much 
smoother overall, which results in a more layered/string-like ap- 
pearance (see Figs. [8] [12] [14] IB.2I and App. As the bow 
shock becomes increasingly circular with larger / (Fig. [12]), the 
instabilities in the slow and fast models appear more, and more 
distinct, with the fluctuations from the latter occurring at much 
larger wavelengths, -0.2 pc. The vortex on the upper right of the 
bow shock is similar to those formed as a result of the K-H in- 
stability and dragged downstream in the adiabatic case. Indeed, 
K-H instabilities are likely to arise in the fast model due to the 
shear between the 73 km s"^ ISM and the 17 km s"^ RSG wind. 



4.1.2. Emission maps 

The strongest radiative emission for both the fast and slow mod- 
els originates just ahead of the contact discontinuity, forming a 
bright ridge along that surface (see Fig. [7] [bottom]). This over- 
all emissivity is a sum of the contributions from 15 different 
coolants (see App. |C] for more details). While some species ra- 
diate from the entire bow shock surface, e.g. AH20(r), others like 
Aco(r) are almost entirely confined to the reverse shock or for- 
ward shock. This has important consequences for the appearance 
of the bow shock shell as shown in Figs. [T3] and [T4] For exam- 
ple, the emission from the forward shock (hotter gas) results in 
a much smoother shell accentuated by small gentle 'hill-like' 
ripples, e.g. Aatomic, whereas emission from the reverse shock 
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Fig. 12. Hydrogen column density (on a logarithmic scale) at 76 000 years for models B [top] and Dl [bottom] seen at inclination 
angles ° [left] to 90 ° [right] in 30 ° intervals@ 



(the shocked RSG wind) produces a layered, more 'cloud-like' 
appearance, e.g. Ah20(H2)- Several coolants, such as gas-grain, 
rotational transitions of CO and H2O, and the heating species 
produce a more 'finger-like', clumpy structure. 

4.2. Luminosity of the bow shock 

When the ISM and RSG winds collide, most of their kinetic en- 
ergy is thermalised. An upper limit for total radiative luminosity 
in this case (assuming n o other energy sources are present) is 
^tot ~ ^wind + ^iSM ( Wilkin et al.|[r9971 where the kinetic lumi- 
nosity of the RSG wind is 

^wind = \m^vI = 2.8 X lO^^ergs s"^ , (11) 

and the ISM luminosity is 

EisM = ^Mwv^ = 1.6 X lO^^^-^ergs s"^ . (12) 

In reality, only a fraction of the kinetic energy will be con- 
verted, therefore £tot is an upper limit for the radiative luminos- 
ity. EisM for the fast model is 5.3x10^^ ergs s"\ and is 1.1x10^^ 
ergss"^ for the slow models assuming an average ISM density 
of nil = 1.5 cm"^. The combined mechanical luminosity £tot 
(see Table. ^ therefore does not exceed ~6 x 10^^ ergss"^ for 
Betelgeuse's parameters. This is a fairly robust upper limit be- 
cause (a) for the most energetic case (model D) the wind kinetic 
luminosity is less than the ISM contribution; (b) uncertainties in 
the observed Bet e lgeuse mass-loss rate are towards lower values, 
e.g. 'Young etaP (Il993ah : and (c) the stellar velocity cannot be 
a factor of two greater since this scenario would produce a very 
bright forward shock at shorter wavelengths (cf. Mira and IRC 
10216 -see Sect. [5]). 

We compare these theoretical values with estimates of the ra- 
diative luminosities from the simulations also given in Table. IH 

^ Movies showing the bow shock rotated through different inclination 
angles are included in the electronic version. 



We used the following procedure to isolate only the contribution 
from the bow shock. For each model, a (1/r^) density profile 
was utilised to subtract the unshocked RSG wind, then only par- 
ticles with both densities and temperatures greater than the am- 
bient ISM were selected. We further restricted the integration to 
within 6 < 115° (x < Rso) the dominant emission region in the 
AKARI and IRAS images (the cut-off' also serves to clearly delin- 
eate an emission region that can be directly compared with fu- 
ture high-resolution observations). The luminosities for the slow 
models do not diff'er significantly (the increasing stellar velocity 
largely compensates for the decreasing ISM density from model 
A to C), thus we only include the values for model B in Table. IH 
The total bow shock luminosity for the fast and slow models is 
approximately 16% and 29 % of the theoretical luminosity, re- 
spectively. Since £tot is the expected luminosity over 4n: sterad 
and we only include emission from x < Rso, these values are a 
lower limit estimate for the fraction of kinetic energy that is ther- 
malised. In Fig. [151 we show both the total radiative luminosities 
given in Table. |4] and the individual ISM and RSG contributions 
of each species. Since the stellar wind properties are the same 
for all the models, it is not surprising that the luminosities pro- 
duced by the RSG wind component are similar for both the fast 
and slow models. In contrast, the ISM contribution in the fast 
and slow models diff'ers significantly: the higher density in the 
slow models leads to greater luminosities for all species except 
for the atomic lines, and the higher shock temperature in the fast 
model results in higher Latomic (Fig. [TSl [bottom]). 

We can estimate the luminosity in the IRAS and AKARI 
wavelength bands from the observed ffuxes: 

Ly = 47rD^FyAy, (13) 

where Ay is the instrument bandwidth at the given wavelength 
and Fy the ffux measured in that wavelength band. This yields 
an underestimate of the bolometric luminosity. The AKARI ffux 
is 10.7 Jy and 6.9 Jy at 65 jim and 90 yum, respectively (lUeta et al.l 
2008, private comm.), and the IRAS ffuxes are significantly 
greater wit h 1 10+20 Jy and 40+1 Jy at 60yum and lOOyum, re- 
spectively (iNoriega-Crespo et al.lll997h . The derived luminosi- 
ties are given in Table. |4l and show that, while the IRAS lumi- 
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Fig. 13. Projected emissivity (ergs cm"^ s"^) in the symmetry plane for model Ah for the various molecular, atomic and dust species. 



nosities exceed the theoretical upper limit, the AKARI values are 
consistent given the uncertainties. Contamination by flux from 
the bar and Betelgeuse itself, which is very luminous in the in- 
frared, may play a role and are likely to account for the difl'erence 
between the AKARI and IRAS values. Higher resolution observa- 
tions, together with a careful subtraction of the bar and central 
source, are essential for resolving this discrepancy. 

The most important coolants at low gas temperatures, 
-300 K, are the gas-grain cooling, rotational modes of CO and 
H2O, and the rotational and vibrational transitions of H2. The 
dust, and C and O fine structure lines are often invoked to ac- 
count for the infrared emission from bow shocks. It is possible 
that including non-LTE chemistry could change the abundance 
of 01, however it appears to contribute little to the far-infrared 
emission. In contrast, the dust and CI, CII fine structure lines ap- 
pear to be the dominant far-infrared emitters in the bow shock 
shell, although it shou ld be noted that the c ooling rate of the lat- 
ter is very uncertain ( Smith & Rosenll2003[ private comm.). The 
combined luminosity derived from the simulations due to grains. 



01, CI, and CII is at least three orders of magnitude lower than 
inferred from the AKARI observations. A possible explanation is 
that there is an additional source of infrared flux, e.g. some of 
Betelgeuse's radiation, or radiation produced by hot gas in the 
bow shock itself, is absorbed and reemitted by the gas and dust 
in the far-infrared. 



5. Discussion 

5.1. Comparison with previous studies 

Our results are largely consistent with the previous hydrody- 
namic studies highlighted i n Sect. [T] The slow model s are 
similar to the 2D mod els of iBrighenti &D'Ercolel(ll995h and 
van Marie et in their flow characteristics and their R- 

T substructure . Howe ver, the bow shocks in the 3D models of 
Wareing et al. (l2007ah are generally smoother, and where insta- 
bilities are present, they occur on much larger scales than in 
our models. The difl'erence is likely due to our including low- 
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Fig. 14. Same as Fig. [13] except for model D. 



temperature cooling (below 10 000 K - the temperature they as- 
sume for the stellar wind) and may also be a resolution effect. 

The smoother bow shock and the prominent vortex in the up- 
per right of the fast model appear to be more consistent with the 
AKARI observations of Betelgeuse's bow shock than the highly 
unstable bow shock in the slow models. However, with such high 
ISM shock temperatures (10^ K), the bow shock should be vis- 
ible at shorter wavelengths, e.g. in the UV. To our knowledge, 
no such emission has been detected for Betelgeuse. For Mira (o 
Ceti), an AGB star, 5 x 10^ K gas was detected serendipitously 
with the GALEX UV satellit e in both the bow shock and a 4 pc 
long tail (iMartin et al.l l2007h . The star is estimated to be moving 
through a low density , ^h~0.03 - 0.8 cm"^ ISM at 130 kms"^ 
(IWareing et al.ll20079: iMartin et al.ll2007l) . The emission mech- 
anism for this UV radiation is uncertain, but is thought to be 
the interaction of hot elec trons (heated in t he forward shock) 
with molecular hydrogen (iMartin et al.l 1200 7). The hot gas in 
the core of the tail for the fast model is certainly reminiscent 
of Mira's narrow tail, and would provide a ready supply of elec- 
trons. Mira's bow shock has also more recently been detected in 



the infrared with Herschel (Mayer e t al .11201 ll) . Another system 
with both UV and confirmed infrared emission from the bow 
shock is IRC 10216, an AGB star movin g through th e ISM at 
~90kms-i (ISahai & Chronopoulosll2010l:lLadial et aDl2010h . 

High-resolution observations, e.g. with Herschel and SOFIA, 
should be able to distinguish between the fast and slow mod- 
els based on the morphology of the bow shock and the pres- 
ence/length scale of the instabilities. 

5.2. The age of Betelgeuse's bow shock 
5.2.1 . Bow shock shell mass 

By comparing the bow shock shell mass derived from the mod- 
els with observational estimates, we can constrain the age of the 
bow shock shell. In Fig. [16] we show the evolution of the bow 
shock shell mass for models B [left] and D [right] . Contributions 
from material in the bow shock were obtained using the same 
procedure as for the luminosity calculations described above, es- 
sentially only shocked particles in the bow shock head (i.e. with 
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Table 4. Bow shock luminosities in ergs s"^ . 



Species 


Model B 


Fast model 


1- ^total 


3.8x10^2 


9.2x10^2 


2. ^grain 


4.1x10^^ 


2.5x10^^ 


3. ^0(63//m) 


2.7x10^5 


2.1x10^6 


4. Lciji 


9.7x10^0 


2.7x10^0 


5. ^H2(r-v) 


8.1x10^1 


1.0x10^2 


6. ^atomic 


2.2x10^2 


7.9x10^2 


7. ^H20(r) 


6.2x10^1 


2.4x10^1 


8. ^H20(H2) 


7.5x10^6 


3.0x10^^ 


9. ^H20(H) 


2.2x10^0 


9.5x10^9 


10. LcO(r) 


5.2x10^^ 


4.8x10^^ 


11. LcO(H2) 


1.3x10^6 


4.8x10^6 


12. Lco(H) 


8.5x10^^ 


1.3x10^^ 


13.LoH 


2.3x10^8 


5.1x10^8 


Theory 




1.3x10^^ 


5.6x10^^ 


Observations 


^AKARI 


6.9x10^^ (65yum) 


2.5x10^^ (90yum) 


^^AKARI 


10.7 Jy (65^m) 


6.9 Jy (90yum)) 


^IRAS 


4.7x10^4 (60yum) 


1.87x10^4 QQQ^j^) 


^^RAS 


110+20 Jy (60yum) 


40+10 Jy(lOOyum)) 



Notes. [1] T. Ueta, private comm. [21 iNoriega-Crespo et al.ll 19971 . 

X < 0, points in Fig. [16]) and x < Rso (solid lines in Fig. [161) were 
included. The former is useful for comparison with the theoret- 
ical estimates, while the latter is utilised when discussing the 
observations. The mass in the bow shock increases almost lin- 
early at early times and only reaches an equilibrium state for 
the bow shock head after at least 60 000 years. The slow models 
take longer to achieve this steady state because the gas dragging 
material downstream has much lower velocity than in the fast 
models. Most of the mass in the bow shock shell is from the 
RSG wind, and as expected, its contribution does not differ sig- 
nificantly between the fast and slow models, except at later times 
when more wind material is able to accumulate at the bow shock 
due to mixing caused by R-T instabilities in the latter. As was 
the case for model Bad, the RSG shell mass for x < is approx- 
imately -0.05 Mo, twice that of the theoretical estimate derived 
using O.SxAtsheii- Owing to the higher ISM density assumed, the 
ISM contribution in model B is larger (half as much as the RSG 
wind) than that in model D (only about a quarter), and as a result 
the total mass in the bow shock head, >O.O8M0, is also greater 
than in model D, >0.048 M©. 

The mass in the bow shock based on the observed 60 jim flux 
is given by 



Atsheii(obs) = 0.042Mo(F6o/135 Jy)(Z)/200pc)2 



(14) 



Noriega-Crespo et aPl 19971) . Again, assuming the lHarper et al.l 
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Fig. 15. Models B and D radiative luminosity contributions for 
the various molecular and atomic species ordered on the x-axis 
as they are listed in Table. [H The total luminosity is shown in 
black [top], with the RSG and ISM components shown in red 
and blue [bottom], respectively. 




Tot, x<R, 
ISM, x<R3Q 
RSG, x<RsQ 
Tot, x<0 
ISM, x<0 
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40 60 20 40 60 80 

Time [1000 years] 

Fig. 16. Evolution of the bow shock shell mass for models B 
[left] and D [right]. For both models, the RSG and ISM con- 
tributions are in red and blue, respectively, and their combined 
mass is plotted in black. The solid lines trace mass for x < Rso, 
while the points trace x <0 gas. 



2008h distance of 197 pc, and the IRAS Mx of 110 Jy in 7.5 



arcmin yields a shell mass of 0.033 M©. (iNoriega-Crespo et aH 



1997" derived a sheU mass of 0.14 M©, but assumed D = 400 
pc for the distance to Betelgeuse.) The corresponding age of 
the bow shock is approximately 30000 years (see Fig. [161 solid 
lines). However, as discussed in Sect. 14.21 this IRAS flux is most 
likely an overestimate, thus the age derived based on that value 
is an upper limit. 

The shell mass from the AKARI fluxes is an order of magni- 
tude lower (-0.0033 M©), which would imply an even younger 
bow shock age. If this is the case, however, the wind would not 
have sufficient time to expand to the stand-off' distance. One pos- 
sible solution is that the observed shell mass is underestimated 
due to uncertainties in the conversion from ffux to mass (e.g. due 
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to dust properties). In our models the shell takes -20000 years 
to reach the correct Rso (see Fig. O, by which time the mass 
in the bow shock is approximately 0.02 M©. This is higher than 
the value obtained from AKARI observations, but may be con- 
sistent within the uncertainties and with the consensus estimates 
for Betelgeuse's wind parameters. However, at this age none of 
our models are close to reaching a steady state. 

Another possibility is that the shell is older; however, this 
requires extreme stellar wind properties. Utilising Eq.|7l we see 
that decreasing the mass-loss rate and increasing the wind veloc- 
ity will both decrease the shell mass. Howe ver, the largest line 
widths observed in the wind are 40kms-i (Hu22ins eiai]ll994l 
Fig. 1), thus the wind velocity, which is usually taken to be half 
of this linewidth, must be approximately 20 km s"^ or less. Thus, 
to decrease the shell mass from a steady state value of 0.05 M© 
to 0.0033 Mq requires that we decrease the mass-loss rate by a 
factor of 15, i. e. M = 2x10 " ^ M© yr ~\ somewhat below the low- 
est estimate of* Youn g et al.l (Il993ah . The latter was derived from 
CO observations, which tends to underestimate the true mass- 
loss rate due to inco mplete C O synthesis in Betelgeuse's stellar 
wind (Noriega-Crespo^et^ [l997). Although this lower mass- 
loss rate and a higher wind velocity cannot be excluded, they are 
outlying values in the range of observations so must be consid- 
ered unlikely. 

5.2.2. The shape of the bow shock 

Based on these arguments, our results suggest that Betelgeuse's 
bow shock is young and may not have reached a steady state yet. 
The circular nature and smoothness of bow shock can be natu- 
rally accounted for if this is the case, and need not be due solely 
to the inclination of the bow shock with respect to the plane 
of the sky; however, Betelgeuse's tangential velocity, yt=25 
kms"\ and radial velocity, yrad=20.7 kms"^ (assuming the dis- 
tance to the star is 197 pc, see Sect.[T]) yield an inclination angle 
of arctan(yt/yrad)~50°. Given the uncertainty in the distance to 
the sta r, this inclination angle is consistent with the S Ueta et af] 
(•2008") value (56°) derived by fitting the shape of the observed 
arc with the analytic Wilkin (1996) solution, even though as dis- 
cussed in Sect. 14.11 the latter is only valid if the bow shock has 
reached a steady state. 

Utilising Betelgeuse's radial and tangential motions, the stel- 
lar velocity is (V^ + y^^^^)^/^~32.5 kms"^ the same as model B. 
Although the smooth appearance of the shell would seem to rule 
out the slow models, if the bow shock is young, the strong in- 
stabilities that characterise those simulations may not have had 
enough time to grow. We can estimate the age of the bow shock 
by comparing the observed shape with what is predicted by the 
simulations in Fig. [TT] [top]. From the AKARI observations, the 
ratio of R(0°)/R(90°) is approximately 0.7, which corresponds 
to an age of <30 000 years. 

According to our simulations, if the bow shock is young, 
the bow shock tail should show strong curvature and should 
not be too distant from the head of the bow shock. Although 
the emission from such a tail is likely to be weak, deep, high- 
resolution observations may be able to detect it. Indeed, there 
appears to be faint emission extending from the bow shock head 
towards the tail in the AKARI observations, e.g. the structure lo- 
cated at (RA,DE C) offsets o f (-2,9) arcmin in the WIDE-S im- 
age of Fig. 1 of lUeta et al.l (^008). A curved tail has already 
been observed for the massive 09.5 supergiant a Camelopardalis 
(S. MandeQ) and likely for the 06 I(n)f runaway star A Cep 

^ See image on APOD at ,http://apod.nasa.gOY/apod/ap061124.ht'inr] 



(iGvaramadze & GualandrisI l201ll) . Although these are much 
hotter stars with faster stellar winds, the physical mechanism 
responsible for producing a curved tail is the same. In light of 
our simulations, such a feature implies either a very small space 
motion for the star relative to ambient ISM (unlikely given the 
strength of the bow shock emission) or that the bow shock is 
young. Thus, modelling the shape of the bow shock (head and 
tail, if present) is a promising avenue for constraining the age 
and the evolutionary stage of these systems. 

5.2.3. Implications 

Estimates for Betelgeuse's age are highly uncertain, ranging 
from 8-13 million years depending on the stellar evolution model 
and distance assumed (iHarper et al.l ^08). The core helium- 
burning phase lasts about 10% of the total lifetime, thus thc = 
0.8 - 1.3 Myr. The star is expected to start and spend a sig- 
nificant fraction of this phase as an RSG. According to cur- 
rent stellar models, the star may or may not have undergone a 
so-called blue loop in the HR diagram, spending a fraction of 
core helium burning as a blue supergia nt, before returning to the 
RSG stage at core heli um exhaustion (iMeynet & Maedeij|2000l: 
iHeger & Langer||2000|). 

A young age of Betelgeuse's bow shock might imply that the 
star has entered the RSG stage only recently. It may be interest- 
ing to relate this idea with the recent finding that the diameter of 
Betelgeuse at 1 1 jim has systematically decre ased by about 15% 
over the p ast 15 years (Towne s et al.l 12009*). While iRavi et all 
( 2010) and lOhnaka et al. ( 2011) debate whether this implies a 
real radius change, or rather a change in the density of certain 
layers in the envelope of the star, it is remarkable that the time 
scale of these changes is of the order of 100 years. This might 
imply that the envelope of Betelgeuse is not in thermal equilib- 
rium, as might be the case when a star is entering the RSG stage. 

On the other hand, the radius of a star entering the RSG 
phase for the first time after core hydrogen exhaustion only in- 
creases, with time scales for the radius change (R/R) down to 
about 1 000 years. This is different, however, when the star enters 
the RSG stage after core helium exhaustion, because the igniting 
helium shell leads to an intermittent episode of shrinkage. Since 
the star might have lost a substantial fraction of its hydrogen-rich 
envelope by that time, its thermal time scale might also have de- 
creased. Thus, while more observational and theoretical efforts 
are required to give this speculation more substance, a consistent 
picture might involve Betelgeuse recently having finished core 
helium burning and returning to the RSG stage from a previous 
blue supergiant excursion. 

If Betelgeuse has only recently entered the RSG stage, it may 
not have had enough time to travel beyond its main sequence or 
blue supergiant wind bubble. Assuming a wind mass-loss rate of 
10"^ M© yr"\ and a wind velocity of ~10^ kms"\ the stand-off 
distance for Betelgeuse's main sequence or blue supergiant bow 
shock shell was around 1 pc. A RSG phase of a -fewxlOOOO 
years would bring the star close to the edge of such a bubble. 
This raises the possibility that the mysterious 'bar' ahead of the 
bow shock could be a remnant of this earlier phase of evolution. 
A theoretical investigation of such a scenario is currently under- 
way. 



6. Conclusions 

We presented the first 3D models of the interaction of 
Betelgeuse's RSG wind with the ISM. We took dust, atomic- 
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, molecular-, and metal-line cooling into account. The mod- 
els cover a range of plausible ISM densities of 0.3 - 1.9 cm"^ 
and stellar velocities of 28 - 73 kms"^ We showed that the 
flow dynamics and morphology of the bow shocks in the mod- 
els difl'ered due to the growth of Rayleigh-Taylor or Kelvin- 
Helmholtz instabilities. The former dominate the slow models, 
resulting in a clumpy substructure, whereas the latter are char- 
acteristic of the fast model and produce a more layered sub- 
structure. In the fast model, gas is shocked to high temperatures. 
If gas as hot as this were to be detected at short wavelengths 
(e.g. UV), this would exclude the slow models as an explanation 
for Betelgeuse's bow shock. High spatial and spectral resolution 
observations (e.g. Herschel, SOFIA and ALMA) particularly of 
the more dominant cooling/emitting species, e.g. rotational lines 
of H2O or CO, could also be used to further constrain the phys- 
ical characteristics of the system. In addition, better determina- 
tion of the molecular to atomic hydrogen fraction in the RSG 
wind, the abundances of various species, and temperature of the 
ISM would also reduce the number of free parameters in the 
model. 

The large fluxes in the infrared compared to the theoret- 
ical limit for the bolometric luminosity suggest that the stel- 
lar flux and/or flux from hotter gas in the bow shock is repro- 
cessed by the dust and reemitted in the far-infrared. We showed 
that, if the bow shock shell mass is low, as is suggested by 
the AKARI fluxes, then Betelgeuse's bow shock is young. The 
smoothness and circular nature of the bow shock would be con- 
sistent with this conclusion. Furthermore, if the bow shock has 
not yet reached a steady state, we are less able to constrain the 
physical parameters of the system, e.g. the ram pressure of the 
ISM. This also raises the intriguing possibility that Betelgeuse 
has only recently become an RSG and that the mysterious 'bar' 
ahead of the bow shock is a remnant of a wind shell created dur- 
ing an earlier phase of evolution. 
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Appendix A: The effect of resolution 

As discussed in Sect. 13.11 it is important to use a large num- 
ber of particles to reduce shock broadening. High resolution is 
also necessary in order to resolve instabilities, particularly for 
the slow models. To test this, the slowest model (A) was rerun 
as model Ah with an average six times higher mass resolution, 
and the results are compared in cross-section in Fig. lA.ll [topi. 
A higher resolution version of model D was not computation- 
ally feasible, so a lower resolution version (model Dl) was run 
instead with an average mass resolution that is two times lower, 
and the results are shown in the lower two panels of Fig. lA.ll 
In both cases, as expected, the higher resolution model shows a 
smoother RSG wind and stronger instabilities developing in the 
bow shock. Model Dl appears to have insufficient resolution to 
fully develop the instabilities seen in model D (e.g. the vortex 
in the upper right of the bow shock), whereas the difference be- 
tween models A and Ah is less dramatic because both models 
show the development of R-T instabilities. Despite these differ- 
ences, the overall agreement between the larger scale fluid prop- 
erties of the high and low resolution simulations is very good, 
and they do not differ by more than a few percent. Thus we 
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Fig.A.l. Cross-sectional density profiles for models Ah [top, 
left], A [top, right], D [bottom, left], and Dl [bottom, right] in 
the symmetry plane. 




Fig.B.l. Temperature [left] and total emissivity [right] in the 
symmetry plane for models Cl65o [top row], Cli600 [second 
row], Clsooo [third row], and CLf [bottom row]. 



conclude that our calculations are sufficiently numerically con- 
verged in terms of the quantities we are interested in measuring. 



Appendix B: The effect of varying Tism and / 

While the thermal pressure of the ISM is insignificant compared 
to its ram pressure in all models, it could aff'ect the bow shock 
thermal and emission properties. The wind thermal pressure is 
insignificant owing to adiabatic expansion, but the hydrogen 
molecular fraction, /, entering the reverse shock may also aff'ect 
the emissivity and cooling of shocked gas. The eff'ects of the ISM 



temperature were tested by using three ISM temperatures for dif- 
ferent runs of model C: Tism = 650 K for model Cl65o, 1 600 K 
for model Cli6oo and 8 000 K for model Clsooo- All models use 
a molecular fraction / = 0.001 (on a scale where is atomic 
and 0.5 is fully molecular) for the wind boundary condition, so a 
comparison model CLf was run with / = 0.18. Cross-sections of 
the gas temperature and total emissivity for these four models are 
shown in Fig. lB.li and can be compared to the default model C 
in Fig. [71 although it should be noted that the logarithmic scales 
are diff'erent in the two figures. Higher ISM temperatures pro- 
duce a smoother forward shock and result in a much fainter tail. 
The hotter gas that accumulates along the symmetry axis in the 
tail is also closer to the reverse shock, filling the void behind the 
RSG wind more efficiently. This is a relatively minor diff'erence 
between the models, however, so the strongest impression we get 
from comparing the top three panels in Fig. IB. II is how similar 
they are in both the structure and emissivity of the bow shock. 

The molecular-to-atomic hydrogen fraction in the RSG wind 
increases as the wind cools and expands away from the star. 
The initial value depends on the temperature structure and 
chemistry near the stellar surface. We assumed a low value 
since Betelgeuse is thought to have a hot, -8 000 K chro- 
mosp here that would dissociate molecules ( Hartmann & Avreti 
However, more recent studies suggest that the tempera- 
ture in this region could be much lower at approximately 4 000 
K (Harper 2010, and references therein), and thus / could be 
higher. This was tested in the / = 0.18 model CLf whose re- 
sults are shown in the bottom panel of Fig. IB. II As expected, 
the cooling is much stronger for larger /, which leads to a more 
unstable bow shock - the apex of the bow shock has moved in- 
wards in model CLf and its shape is far from a smooth curve 
for both the forward and reverse shock. The regular appearance 
and relative stabi lity of the observed bow shock on similar scales 
(lUeta et al .1120081) suggests that the hydrogen is indeed mostly in 
atomic form. 



Appendix C: Emission from the bow shocl< shell 

The emissivities for each of the cooling and heating species in- 
cluded in the simulations are shown in cross-section in Fig. IB. 21 
for models Ah [top] and D [bottom]. These cross-sectional pro- 
files show the location and morphology of the emitting regions 
for each species more clearly than the projections shown in 
Figs.[T3]and[T4l Here we also show results for each coolant indi- 
vidually. As expected, the processes that involve molecular hy- 
drogen, e.g. AH2(r-v), tend to emit strongly from the gas in the 
reverse shock. The 01 fine structure and CO rotational cooling 
are also more concentrated in the denser, cooler regions of the 
reverse shock. The other species tend to emit more evenly from 
the entire bow shock, with the atomic cooling greatest in the for- 
ward shock. In the slow models, the atomic cooling is largely 
absent from the cool, dense R-T 'fingers'. In the fast model, 
several species that do emit from both the forward and reverse 
shocks tend to show a strong ridge of emission just ahead of the 
contact discontinuity, e.g. the grain, water, atomic, and carbon 
fine- structure cooling. The lower ISM density in the fast model 
results in much weaker emission from the bow shock tail than in 
the slow models, e.g. Ah20(H2) and Aco(H2). 
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Fig. C.l. Snapshots of the movies for model B [left] and D [right] simulations. First, the evolution of the gas column density is 
followed, then the final structure is rotated 360°, demonstrating the changing morphology of the bow shock with diff'erent inclination 
angles. 



